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Abstract 

We report the results of large scale numerical studies of the dynamic and 
static properties of the random phase model of the randomly pinned flux 
lines confined in a plane. From the onset of nonlinear IV characteristics, we 
identify a vortex glass transition at the theoretically anticipated temperature. 
However, the vortex glass phase itself is much less glassy than expected. No 
signature of the glass transition has been detected in the static quantities 
measured. 
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Flux pinning is crucial to the performance of high-T c superconductors in a strong mag- 
netic field It has been conjectured [FJ,[3| that flux lines in a superconductor with random, 
point-like impurities may form a glass, the vortex glass, in which the flux line configurations 
are frozen by the defects at low temperatures. It is argued ||rj]] that if a glass phase ex- 
ists, then the sluggish glassy dynamics of the flux lines suppresses dissipation, making the 
system a true superconductor. Unfortunately, evidence supporting the vortex glass hypoth- 
esis is still inconclusive: The phase transition seen experimentally by Koch et al || and by 
Gammel et al || may be due to pinning by correlated defects such as twin planes and/or 
screw dislocations. Existing numerical simulations which support the existence of the vortex 
glass are thus far restricted to models with zero external magnetic field and/or an infinite 
London penetration length ||. However, the glass phase seems to be unstable to thermal 
fluctuations in more realistic models with a finite London penetration length || . 

Without solid evidence from experimental and numerical studies, major support for the 
vortex glass hypothesis comes from a number of analytic studies [p|J^ JT0|JTT|1 , which are how- 



ever restricted to "elastic systems" where topological defects are forbidden to nucleate. While 
the validity of the elastic approximation for the bulk system in 2 + 1 dimensions is a subject 
of intense current research, the elastic model is well justified to describe flux lines confined 
in a plane (1 + 1 dimensions) j2l Jl"2yi3H or vortex lines in planar Josephson junctions [ 14| . 



Thus far, the 1 + 1 dimensional flux array is the only system for which a vortex glass phase 
has been demonstrated analytically; as such, it is one of the very few pieces of "solid" sup- 
port on which the vortex glass hypothesis for bulk superconductors is based. Because of 
this significance, the static and dynamic properties of the 1 + 1 dimensional system have 



been investigated analytically by many groups in the past several years P,|T0|-|17|, ? ,|19H22 
These studies led to a number of rather different predictions, none of which has been tested 
experimentally or numerically. 

In this paper, we report the results of detailed numerical studies of the random phase 
model P JT^JT^J2B| of the planar flux array using the Connection Machine CM5. We ob- 



serve an apparent vortex glass transition at a finite temperature, T g , below which the IV 



characteristics are nonlinear and described by power laws, as anticipated by Refs. |L2| and 



T6| . However, the universal constant found is significantly smaller than that predicted in 



Ref. (T^|. Static quantities such as the equal-time correlation function are also measured. 
They exhibit no detectable signs of the glass transition. 

We consider an array of flux lines confined to the (x, 2;)-plane, and directed in the z- 
direction by an applied magnetic field H = H$z [f23|l - Line-line repulsion is modeled by 
linear elasticity, with an elastic constant k = (dp/dHo)~ 1 , where p is the average line density. 
Random point-like pinning centers in the plane are described by an uncorrelated Gaussian 
random potential with a variance Ao. A direct simulation of the flux array is deemed difficult 
because an interline spacing of ~ 10 grid points and a system size of ^> 10p _1 is needed 
to probe the asymptotic behavior. Instead we will study the random phase model, 

H = J dV{|(Vu) 2 -Acos[27ra(r)-/3(r)]j, (1) 



which is generally believed p. |T2| . |23| to describe the statistical mechanics of the random flux 
array at length scales large compare to the inter-line spacing p~ l . In Eq. ([!]), r = {x,z}, 
u(r) is a displacement-like field, the cosine term captures the discreteness of the flux lines 
crucial to pinning, j3(r) is a random phase uniformly distributed in the interval [0, 2tt], and 
A = 2y/Aop characterizes the strength of the pinning potential. A derivation of this model 
starting from a system of directed flux lines can be found in Ref. [ [23|j . 

Renormalization group (RG) studies of Eq. ([I]) p^j25|| found a glass transition at T g = 
k,/tt for small A. For T > T g , randomness (the cosine term) is irrelevant, and the system is 
asymptotically described by Eq. ([[]) with A = 0. But for T < T g , the system is described 
by a line of fixed points. A number of properties of the glass phase close to T g have been 
obtained from the RG studies. For instance, the static correlation function is predicted to 
diverge with a logarithmic anomaly ||25|| , 



C(r) = ((«(r) - u(0))2) = AC (r) + Br 2 [\n \r\] 2 } (2) 
where r = 1 — T/T g , Co{r) = T/(ttk) In \r\ is the correlation function with A = 0, A is a 



nonuniversal constant, and B = 2/ir 2 + 0{r) is universal | 22| . Overbar and (...) denote 
disorder and thermal averages respectively. 

Dynamic RG methods were used to extract the glassy dynamics of the model ([1]) close 
to T g |T6| . |26| . In particular, under a small driving force F, the Langevin dynamics 

dtU = J2i + F + v , (3) 

ou 

where rj(r,t) is a white noise with (rj(r, t)r](0, 0)) = 2T8 2 (r)5(t), is expected to yield a 



nontrivial response, v (F) = (d t u), with 

v(F) oc |if F, for T > T g , (4) 

v{F) oc |r| c F a , a = 1 + Cr for T < T g , (5) 

with the universal constant ( ~ 1.78+0(r). F describes the effect of a Lorentz force resulting 
from an electric current perpendicular to the (x, z)-plane. v(F), the "IV characteristics", is 
proportional to the EMF generated from the in-plane motion of the flux lines when pushed 
by the Lorentz force |13|]. A zero DC resistivity (limF^ov(F)/F — > 0) and hence true 
superconductivity is expected in the vortex glass phase according to Eq. (H). 

Within the RG framework, these predictions are valid in the asymptotic limits of the 
system size L — > oo and F — > 0. Finite size effects are derived using the length dependence 
of the dimensionless coupling constant g(L), 

_ g(l)L^ 
9[ ' l + Cg{l)[L^ -l]/(2r)' [ > 

which characterizes the effective strength of disorder at scale L, and is obtained from the 
RG recursion relation p2| , p5| . In Eq. (||), g(l) = n\ 2 /(2T 2 ) is the bare coupling constant, 



and the nonuniversal coefficient C ~ 1 controls the crossover length for the model ([[]). 
Finite size effects in the correlation function shall be estimated by replacing r in Eq. (0) by 
(2/C)g(L) 12^. 

When F^O, the equilibration time is t eq ~ 1/v (F), since for v(F)t 3> 1, (u) ~ v(F)t 3> 
1, which averages out the pinning effect of the cosine [|T6 1 . But during the time t eq , only 



regions of size L cS « (z = 2a [[U|) can be equilibrated. Hence the effective system 
size at a finite drive is L e R « p- 1 / 2 \ While the effect of a large F involves nonequilibrium 
dynamics [16[] and is complicated, at the level of linear response, we shall again estimate the 
effect of a finite F by replacing the coefficient |r|^ in Eq. (|5|) by [(2/C)g(F~ 1 ^ 2 )] 1 ^. 

Although all existing analytic studies agree that T g = k/tt, some results of the RG 



analysis disagree with results by other methods. For instance, Toner argued that the 
form of the IV, for T < T g , should be 



v(F) ~ F e- c v /l °s( F o/^). ( 7 ) 



Also, Balents and Kardar |19j obtained, at T = T g , a much more strongly divergent corre- 
lation function 

C(r) ~ |r|, (8) 

using a Bethe ansatz solution of the random flux array. On the other hand, a similar 
study by Tsvelik |?D| yielded C(r) oc In \ r\ when a different order of thermodynamic limits 
was taken. One criticism of the RG analysis is that the possibility of "replica symmetry 
breaking" , often associated with the existence of meta-stable states, has not been considered. 
A recent variational ansatz with hierarchical replica symmetry breaking scheme [H],pT[| yields 
a correlation function even less divergent than Eq. @, 

C(r) = -£ In | r l ( 9 ) 

UK 

for all T < T g and \r\ 3> 1. However, an earlier attempt [Hj yielded Eq. (]8|). 



In view of the myriads of predictions for this system, we perform a large scale numerical 
study. As we shall be interested in the dynamics, we directly simulate the Langevin equation 
(0) of the model (|I|) on an L x L square lattice with periodic boundary conditions in both the 
x and z directions. To do this numerical integration, we approximate the time derivative in 
Eq. by a finite difference, and use a two step, single noise Runge-Kutta integration algo- 
rithm [p8H 3(|- This algorithm yields results accurate to 0(dt 2 ), where dt is the the discrete 



time step. We ensure that the systematic dt 2 errors are as predicted and under control by 
doing runs at various dfs (with 1024 realizations each) and various temperatures [P0|] . The 
simulations reported below are for runs with an optimized dt = 0.2, and sizes L = 64, 128. 
The parameters are fixed at k — 1.0 and A = 0.15. The IV curves are obtained by varying 
F in the range 0.005 to 0.1, and the glass transition probed by varying T from 0.5 down 
to 0.25. The number of realizations over which we averaged ranged from 32 (L = 128) to 
1024 (L = 64) and is specified where relevant. Note that the range of parameters are chosen 
such that g(l) > g(L) > g(oc) according to Eq. (|). Thus the effective strength of disorder 
should appear stronger than its asymptotic value due to the finite size effect. However, this 
effect is not expected to be large, since the crossover length, L c , defined by equating the two 
terms of the denominator of Eq. (|6|), is only L c = e^m pa 15 at its largest (when r = 0) ||30|| . 

We start from random initial conditions in u. Depending on the sizes and temperatures 
of the systems, we typically use 80, 000 to 200, 000 time steps to equilibrate, and then 60, 000 
to 150, 000 steps to perform measurements. "Thermal averages" are carried out by time- 
averaging over the measurement period for each realization. In all cases, equilibration is 
assured by monitoring the susceptibility, \ — (U d 2 r d x u) 2 ) / (TL X L Z ), whose mean is 1/k, 
at all T, within statistical and dt 2 errors [T2^,^]. As extra precaution, we did some very long 



runs (10 6 steps, 160 realizations) for the L = 64 system at T = 0.285. All measurements 
there agree with the shorter runs. 

We begin by examining the IV curves v (F). Figure 1 shows the IV curves for a 64 x 64 
lattice for various temperatures on a log-log plot. (Note that the vertical axis is divided by 
F to emphasize the nonlinear part of v(F).) According to Eqs. (|j) and (|5]), we expect linear 
IV curves only above T g . From Fig. 1, we see that the data for T = 0.4 satisfy a linear 
relation [|3T], while for T = 0.3 and below, they are nonlinear and well described by a power 
law as in Eq. (^). The straightness of the data strongly biases against the prediction (^). 
The L = 128 lattice yields similar results. Fitting Eq. ([5]) to the IV curves, we can obtain 
the critical temperature, T g , and the coefficient (. From Fig. 2 we see that the exponent a 
indeed varies linearly with the reduced temperature r up to at least r ~ 0.21 (or T = 0.25). 



Chi-squared fits to our IV data for T = 0.25, 0.27, 0.285 and 0.3 (L = 64) and T = 0.25, 0.275 
and 0.3 (L = 128), give T g = 0.32 ± 0.01 and C = 0.15 ± 0.02 for the L = 64 system, and 
T g = 0.34 ±0.01 with £ = 0.15 ±0.02 for the L = 128 system. The close agreement between 
the two systems indicates that our result is not limited by the system size L. Our value for 
T g agrees with the predicted value of 1/tt « 0.318, but the predicted value for the universal 
coefficient ( is an order of magnitude too large. For comparison, we indicate the predicted 
slope of the IV curve at T = 0.25 as the dashed line in Fig. 1, and the predicted a(r) as the 
dashed line in Fig. 2. Clearly, neither can be accommodated by our numerical results. 

One may attempt to attribute the discrepancy between the numerics and the RG pre- 
diction to crossover effects, since even the smallest value of F used corresponds to only 
L c g ^15. However, since v (F) oc [g{F- 1 ' 2 )fF a and g{F- 1 / 2 ) increases with F for the 
range of parameters used, one expects the effective values of a observed to be somewhat 
larger than the expected asymptotic values [p0|1 . It is therefore not possible to explain the 



observed discrepancy by such crossover effects. In fact, it is possible (although unlikely) that 
the observed nonlinear IV's are the pre-asymptotics of some linear IV. A study at smaller 
F's is desirable but difficult numerically, since a much longer measurement period will be 
needed to obtain v(F) accurately for each F. 

Next we consider the static correlation function, C(r), for which we have much better 
statistics and smaller finite size effects since F — 0. In comparing our data for C(r) with the 
predicted forms, Eqs. (|]), and we do not compare directly with In |r| since that is 
the form for the quadratic theory in the continuum. Instead we compare with the correlation 
function Co(r) of the quadratic lattice model calculated exactly for dt = 0.2. Figure 3 shows 
the measured values of C(r) at T = 0.3, 0.285 and 0.25 for the L = 64 system each averaged 
over 1024 realizations. (The error bars are much smaller than the points and are not shown.) 
The data are well described by C(r) = ACo(r) with A ~ 1.09 for all T, and we find A — > 1 
as dt — ► for all r [30|| . Our results are significantly different from Eq. (0), whose prediction 
for C{r) at T = 0.25 is plotted as the dashed line in Fig. 3 and should be easily detectable, 
unless the universal coefficient B is much (two orders of magnitude) smaller than predicted. 



Again, it will be difficult to attribute the discrepancies to finite size effects which would have 
made C(r) appear slightly larger since g(l) > g(oo) f50fl . Our data are also in disagreement 
with Eqs. (§) and @. However, the crossover length to Eq. is expected to be long 
[0(e 1//r )] for small r [II]. Here, simulations at even lower temperatures are needed to make 
an unequivocal test. 

Finally we examined the variance of the susceptibility x, which is expected to be nonzero 



and universal in the glass phase ||22|| . However, our numerical results yield a variance that is 
zero within statistical error above and below T g [|30|| . It is again in clear disagreement with 
the RG prediction and in line with the very small effect of disorder observed in all of the 
other measurements. 

To summarize, our numerical simulation of the random phase model of the randomly 
pinned planar flux array leads to some rather surprising results: Although the onset of 
nonlinear IV characteristics is observed at the expected glass transition temperature, the 
low temperature phase is found to be much less glassy than that predicted by the RG 
analysis. Both the disorder-averaged correlation function and the normalized variance of 
susceptibility variations are, within statistical errors, indistinguishable from those of the pure 
system. Our numerical results clearly rule out the RG predictions — either the universal 
constants ( and B are respectively at least one and two orders of magnitude smaller, or the 
crossover scale is at least two orders of magnitude longer. Other possible scenaria include: 
(1) We are prevented from reaching the asymptotics of some unknown fixed point (such as 
those proposed in Refs. |TT| , p!9| -pT|) by a very strong crossover effect. (2) The phase transition 
observed is purely dynamic and does not show up in the statics. The source of discrepancies 
between theory and simulation is very puzzling and warrants more detailed theoretical and 
numerical investigation. 
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FIGURES 



FIG. 1. The nonlinear part of the IV curves for L = 64 at T = 0.4, 0.3, 0.25. The solid lines are 
best straight line fits. The dashed line is Eq. (5) at T = 0.25. 



FIG. 2. The exponent a for L = 64 as a function of the reduced temperature. The dashed line 
is a = 1 + 1.78t. 



FIG. 3. The correlation function for L = 64 at T = 0.3, 0.285, 0.25 with 1024 realizations each. 
The solid lines are C(r) = 1.09C (r). The dashed line is Eq. (2) at T = 0.25 and A=l. 
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